AddBiomechanics: Automating model scaling, inverse kinematics, and inverse dynamics from human motion data through sequential optimization

Creating large-scale public datasets of human motion biomechanics could unlock data-driven breakthroughs in our understanding of human motion, neuromuscular diseases, and assistive devices. However, the manual effort currently required to process motion capture data and quantify the kinematics and dynamics of movement is costly and limits the collection and sharing of large-scale biomechanical datasets. We present a method, called AddBiomechanics, to automate and standardize the quantification of human movement dynamics from motion capture data. We use linear methods followed by a non-convex bilevel optimization to scale the body segments of a musculoskeletal model, register the locations of optical markers placed on an experimental subject to the markers on a musculoskeletal model, and compute body segment kinematics given trajectories of experimental markers during a motion. We then apply a linear method followed by another non-convex optimization to find body segment masses and fine tune kinematics to minimize residual forces given corresponding trajectories of ground reaction forces. The optimization approach requires approximately 3-5 minutes to determine a subject’s skeleton dimensions and motion kinematics, and less than 30 minutes of computation to also determine dynamically consistent skeleton inertia properties and fine-tuned kinematics and kinetics, compared with about one day of manual work for a human expert. We used AddBiomechanics to automatically reconstruct joint angle and torque trajectories from previously published multi-activity datasets, achieving close correspondence to expert-calculated values, marker root-mean-square errors less than 2 cm, and residual force magnitudes smaller than 2% of peak external force. Finally, we confirmed that AddBiomechanics accurately reproduced joint kinematics and kinetics from synthetic walking data with low marker error and residual loads. We have published the algorithm as an open source cloud service at AddBiomechanics.org, which is available at no cost and asks that users agree to share processed and de-identified data with the community. As of this writing, hundreds of researchers have used the prototype tool to process and share about ten thousand motion files from about one thousand experimental subjects. Reducing the barriers to processing and sharing high-quality human motion biomechanics data will enable more people to use state-of-the-art biomechanical analysis, do so at lower cost, and share larger and more accurate datasets.

In our implementation, the joint jerks ∂ 3 ∂t 3 q are computed using finite differences.

Constructing the b vector for angular dynamics fitting
The vector b ∈ R 6T represents the current COM trajectory through space (first half, 3T entries), and the current root (e.g., pelvis) angular trajectory (second half, 3T entries), at the initial conditions ζ.
To compute z t , note that we are holding the mass of the subject constant during this optimization step (m), so we can simply integrate the effects of known external forces (gravity, GRFs) on the known mass of the subject over time.zt = To compute θ t , we linearly integrate the "residual free angular acceleration" over time.We define the "residual free angular acceleration" with joint state q t , qt , qt to be the necessary root angular acceleration θt such that no angular residual force is present.By convention θt is always the first three entries of qt , which we write qt, [1:3] .
We can compute θt given q t , qt , qt by first solving for joint torques τ t using inverse dynamics.Then, set the first three entries of τ t to 0, and solve forward dynamics using q t , qt , τ t .The first three entries of the resulting qt are the "residual free angular acceleration" θt .
Given the "residual free angular acceleration" θt , we can compute θ t by linearly integrating the "residual free angular acceleration" over time:

Accounting for force plate registration errors
The location of experimental force plates are often registered incorrectly in the motion capture volume, which leads to solutions with center of mass trajectories offset slightly from the experimental marker data.To align our solution with the observed experimental data, we shift both the location of the force plates and the marker trajectories by the offset between the center of mass trajectory we find and the experimentally estimated center of mass trajectory, which typically is only a fraction of a centimeter.It is also common to have ground reaction force data recorded from force plates that are very slightly (e.g., less than 0.5 degrees) off of perfectly vertical.This makes finding a physically-consistent solution challenging, because if we assume that the force plates are perfectly vertical in our optimization problem, the total ground reaction force vector will be very slightly off of perfectly vertical in the ground reference frame.This can lead to a substantial horizontal acceleration bias on very long trajectories, since the horizontal ground reaction forces will not be exactly anti-parallel to gravity.
Since the force plate rotation errors are typically very small, we can preserve the linearity of our system by using a first-order Taylor expansion to approximate the rotation of measured forces by very small angles.In our implementation, we append a rotational correction term α i ∈ R 3 to ζ for each force plate i, where n is the number of force plates.
June 15, 2023 2/4 For each α i , we append to A a 3T × 3 columnar block: where f i t is the ground reaction force vector associated with force plate i.The [.] operator makes a skew-symmetric matrix out of a vector in R 3 , so that a × b = [a]b.
Finally, we regularize to angles α i to discourage large force plate rotations.With this extension to our linear fitting problem, we can recover force plate rotations to a very small fraction of a degree.

Quantifying the impact of using the bilevel constraint in kinematics fitting
We ran an ablation study to quantify the impact of the formulating our kinematics fitting problem as a bilevel optimization problem.Here, we compare two conditions: optimizing a bilevel objective and optimizing a simple monolevel objective.The objective for the bilevel optimization problem is as follows: max s,p To convert the bilevel objective to a monolevel objective, we move joint angle optimization, max qt , to the outer objective along with the body scale and marker registration optimization: We find that using a bilevel approach leads to much faster convergence, and we can get lower marker RMSE at the cost of slightly lower anthropometric prior probability.The bilevel optimizer takes slightly more wall-clock time per iteration, but because it is able to reach high quality marker RMSE in many fewer iterations, it is able to save wall-clock time overall.We ran both optimization functions for several different fixed numbers of iterations on a single walking trial on a commodity server; these results are summarized in Table A1.
max s,p,qt P x( xt |q t , s, p) • P s (s) • P p (p| p)